install.packages("Seurat")
library(Seurat)  
library(dplyr)
library(ggplot2)
library(patchwork)

# 设置工作目录  
setwd("C:\\Users\\29228\\Desktop\\2024春学期各项汇报\\R语言与生物信息学") 
# 解压并加载数据
data_dir <- "pbmc3k_filtered_gene_bc_matrices.tar.gz"
# 设置数据目录
data_dir <- "pbmc3k_filtered_gene_bc_matrices/filtered_gene_bc_matrices/hg19"

# 读取10X数据
data <- Read10X(data.dir = data_dir)

# 创建Seurat对象，初步过滤低质量细胞和基因
pbmc <- CreateSeuratObject(counts = data, project = "pbmc3k", min.cells = 3, min.features = 200)
# 计算线粒体基因比例
pbmc[["percent.mt"]] <- PercentageFeatureSet(pbmc, pattern = "^MT-")

# 根据质控标准进行过滤
pbmc <- subset(pbmc, subset = nFeature_RNA > 200 & nFeature_RNA < 2500 & percent.mt < 5)
# 数据标准化
pbmc <- NormalizeData(pbmc)

# 选择高变基因
pbmc <- FindVariableFeatures(pbmc, selection.method = "vst", nfeatures = 2000)
# 对标准化数据进行缩放
pbmc <- ScaleData(pbmc, features = rownames(pbmc))

# PCA降维
pbmc <- RunPCA(pbmc, features = VariableFeatures(object = pbmc))

# 查找邻居并进行聚类
pbmc <- FindNeighbors(pbmc, dims = 1:10)
pbmc <- FindClusters(pbmc, resolution = 0.5)
# 运行UMAP降维
pbmc <- RunUMAP(pbmc, dims = 1:10)

# 绘制UMAP图
DimPlot(pbmc, reduction = "umap")
umap_plot <- DimPlot(pbmc, reduction = "umap")
ggsave(filename = "umap_plot.png", plot = umap_plot,width = 10, height = 8, dpi = 300)
# 鉴定marker基因
markers <- FindAllMarkers(pbmc, only.pos = TRUE, min.pct = 0.25, logfc.threshold = 0.25)

# 每个cluster的Top2基因
top2_markers <- markers %>% group_by(cluster) %>% top_n(2, avg_log2FC)

# 绘制气泡图
bubble_plot <- DotPlot(pbmc, features = top2_markers$gene) + RotatedAxis()
print(bubble_plot)
ggsave(filename = "bubble_plot.png", plot = bubble_plot,width = 10, height = 8, dpi = 300)
# 根据鉴定的marker基因挑选几个基因进行可视化
feature_genes <- c("CD3D", "MS4A1", "LYZ", "CD8A", "GNLY")  # 示例基因，可以根据实际数据调整

# 绘制Feature Plot
FeaturePlot(pbmc, features = feature_genes)
feature_plot <- FeaturePlot(pbmc, features = feature_genes)
ggsave(filename = "feature_plot.png", plot = feature_plot,width = 10, height = 8, dpi = 300)

# 绘制小提琴图
VlnPlot(pbmc, features = feature_genes)
vln_plot <- VlnPlot(pbmc, features = feature_genes)
ggsave(filename = "vln_plot.png", plot = vln_plot,width = 10, height = 8, dpi = 300)

